Variational methods with coupled Gaussian functions for Bose-Einstein condensates 

with long-range interactions. II. Applications 
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Bose-Einstein condensates with an attractive 1/r interaction and with dipole-dipole interaction are 
investigated in the framework of the Gaussian variational ansatz introduced by S. Rau, J. Main, and 
G. Wunner [Phys. Rev. A, submitted]. We demonstrate that the method of coupled Gaussian wave 
packets is a full-fledged alternative to direct numerical solutions of the Gross-Pitaevskii equation, 
or even superior in that coupled Gaussians are capable of producing both, stable and unstable 
states of the Gross-Pitaevskii equation, and thus of giving access to yet unexplored regions of the 
space of solutions of the Gross-Pitaevskii equation. As an alternative to numerical solutions of 
the Bogoliubov-de Gennes equations, the stability of the stationary condensate wave functions is 
investigated by analyzing the stability properties of the dynamical equations of motion for the 
Gaussian variational parameters in the local vicinity of the stationary fixed points. For blood-cell- 
shaped dipolar condensates it is shown that on the route to collapse the condensate passes through 
a pitchfork bifurcation, where the ground state itself turns unstable, before it finally vanishes in a 
tangent bifurcation. 

PACS numbers: 67.85.-d, 03.75.Hh, 05.30.Jp, 05.45.-a 



I. INTRODUCTION 

In the previous paper [JJ variational methods with cou- 
pled Gaussian functions for Bose-Einstein condensates 
with long-range interactions have been developed. The 
purpose of this paper is to demonstrate the power of this 
approach by applying the variational methods to two dif- 
ferent types of condensates, viz. a monopolar condensate 
with an attractive gravity-like 1/r interaction and a dipo- 
lar condensate. 

Monopolar condensates with an attractive 1/r interac- 
tion, originally proposed by O'Dell et al. [2J, have unique 
stability properties. For a wide range of the scattering 
length the condensate is stable without an external trap. 
Additionally, the gravity-like interaction of a monopolar 
condensate may be an opportunity to investigate usually 
large-scale physical properties like, e.g., boson stars [3] at 
a laboratory scale. The "monopolar" interaction of two 
neutral atoms with positions r and r' is induced by the 
presence of external electromagnetic radiation. O'Dell et 
al. [2J suggest six triads of orthogonal laser beams to in- 
duce the interatomic potential Wi r (r,r') = —u/\i — r'\, 
where u depends on the intensity and wave vector of 
the radiation, and on the isotropic polarizability of the 
atoms. Although this system has not yet been realized 
experimentally it has already served as a model to com- 
pare results obtained analytically and with exact numer- 
ical techniques [IHS]- 

By contrast, Bose-Einstein condensates with a long- 
range dipole-dipole interaction W\ r (r,r') ~ (1 — 
3 cos 2 0)| J*— r'| -3 have been obtained experimentally with 
52 Cr atoms in 2005 by Griesmaier et al. [7J [5], and more 
recently in 2008 by Beaufils et al. [5]. The collapse of 
the condensate also has been subject to extensive exper- 
imental studies [10]. Theoretical investigations have so 
far mostly been based on lattice simulations [TTHT5] or 



on a simple variational approach with a Gaussian type 
orbital [H] • For a review on dipolar condensates see [TS] . 
In this paper we extend and elaborate in more detail 
preliminary work presented in [16j . In Sec. O the re- 
sults for the monopolar condensates are presented and 
discussed. It is shown that only three to five coupled 
Gaussians are sufficient to achieve convergence of the 
mean-field energy, the chemical potential, and the lowest 
eigenvalues of the stability matrix. In Sec. |III| the method 
of coupled Gaussians is applied to dipolar BEC to clar- 
ify the theoretical nature of the collapse mechanism of 
blood-cell shaped condensates. On the route to collapse 
the condensates passes through a pitchfork bifurcation, 
where the ground state itself turns unstable, before it fi- 
nally vanishes in a tangent bifurcation. Conclusions are 
drawn in Sec. IIVI 



II. MONOPOLAR CONDENSATES 

The time-independent Gross-Pitaevskii equation 
(GPE) for a self-trapped condensate with a short-range 
contact interaction with scattering length a and a 
long-range monopolar interaction reads 
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where the natural "atomic" units introduced in Rcf. [4 
were used. Lengths are measured in units of a "Bohr 
radius" a u = h 2 /(mu), energies in units of a "Rydberg 
energy" E u = /i 2 /(2ma 2 ), and times in units of t u — 
h/E u , where u determines the strength of the atom-atom- 
interaction [2], and m is the mass of one boson. The 
number of bosons N can be eliminated from Eq. (nl by 



using scaling properties of the system [H [5] . We define 
r = r/N, tp — iV 3 / 2 i/J, which leaves the norm of the wave 
function invariant, a sc — N 2 a/a u , fl — /i/JV 2 , omit the 
tilde in what follows, substitute fi — > i(d/d£), and finally 
obtain the time-dependent GPE in scaled "atomic" units 
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It is known that for Bose-Einstein condensates with at- 
tractive 1/r interaction two real radially symmetric so- 
lutions, the ground state and a collectively excited state, 
exist. By varying the scattering length a sc , they are cre- 
ated in a tangent bifurcation at a critical value of a sc 
@HB]- Below the tangent bifurcation no stationary so- 
lutions exist. Approaching the bifurcation point from 
higher scattering lengths the condensate collapses. A 
similar behavior is found for dipolar condensates [14] and 
condensates without long-range interactions [17] . 

For the monopolar GPE numerically exact solutions 
exist. The numerical procedure for the direct integration 
of the GPE is introduced in Refs. [UJB]. ^ is our pur- 
pose to investigate the two states connected via the tan- 
gent bifurcation with the coupled Gaussian wave packet 
method. We use the ansatz 



monopolar potential, and the matrix elements read 
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for condensates with spherical symmetry. Following the 
procedure outlined in pQ the dynamical equations for the 
variational parameters read 



7* = 6ia k - v k 



a k = -A{a k Y - -V k , 
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A. Stationary solutions 

Stationary solutions of the GPE can be computed via 
a nonlinear root search of the dynamical equations Q as 
outlined in pQ. The Gaussian parameters obtained from 
the root search are then used to calculate the mean field 
energy 
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for k = 1, . . . , N. The quantities Vq and V 2 fe are obtained 
from the (2JV x 2JV)-dimensional linear set of equations 
(k,l = l,...,N) 
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and the chemical potential 
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where V^ff = 14c + Kn is the sum of the contact and 



For the variational ansatz with only one single Gaussian 
function (JV = 1) the results can be obtained analytically 



(a) -0.02 



(b) 




FIG. 1. (a) Mean field energy i? m f for self-trapped conden- 
sates with attractive 1/r interaction as a function of the scat- 
tering length obtained by using up to five coupled Gaussian 
wave packets (N — 1, . . . , 5, respectively) in comparison with 
the result of the accurate numerical solution of the stationary 
Gross-Pitaevskii equation (exact), (b) Same as (a) but for 
the chemical potential fi. 
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for the chemical potential. 

The results of the variational computations with N = 1 
to iV = 5 coupled Gaussian functions are presented 
in Fig. [Ha) for the mean field energy and in (b) for 
the chemical potential. For comparison the results of 
the exact solution obtained by direct numerical integra- 
tion of the GPE are shown as red triangles. Although 
the variational solution with a single Gaussian signifi- 
cantly differs from the exact calculation, the qualitative 
behavior is the same: two solutions emerge in a tan- 
gent bifurcation at a critical scattering length, which 



is afc v = -3tt/8 = -1.178097 and a£> n = -1.025147 
for the variational and exact calculation, respectively [4] . 
Note that the excited state with higher mean field energy 
has a lower chemical potential than the ground state. 

The main purpose of Fig. [I] is to demonstrate how the 
coupling of only N — 2 to N = 5 Gaussian functions 
drastically improves the results obtained with a simple 
Gaussian type orbital. The coupling of only two Gaus- 
sian functions already leads to a significant improvement 
for both, the mean field energy value and the chemical 
potential. For three or more coupled Gaussians, the re- 
sults in Fig. [I] can no longer be distinguished from the 
numerically exact solution. The enlargements in Fig. [T] 
illustrate the rapid convergence of the results with in- 
creasing number N of coupled Gaussians. 

Detailed values for the mean field energy and the chem- 
ical potential of the ground state and the excited state 
at scattering length a sc = —1 close to the tangent bi- 
furcation are presented in Table [I] The coupling of three 
Gaussian functions already yields an accuracy of more 
than four digits for the mean field energy of the ground 
state. Using up to five Gaussians, we can safely assume 
the variational result as to be fully converged to the ex- 
act numerical computation marked as N — oo. The rapid 
convergence of the critical scattering length a^. with in- 
creasing number of coupled Gaussians is also illustrated 
in Table |U 

In previous work [3] it has been noted that the peak 
amplitude of the exact wave function at the center r = 
is poorly reproduced by a Gaussian type orbital. Here, 
we illustrate that superpositions of Gaussian functions 
can accurately approximate the exact wave function and 
thus also provide the correct peak amplitude. 

In Fig. [2] we choose the scattering length a sc = — 1 close 
to the exact numerical bifurcation at a^ = —1.02515 to 
compare the wave function of the variational and the nu- 
merically exact ground and excited state. As a second 
example, we choose a sc = —0.6 in Fig. [31 which lies far- 
ther away from the critical scattering length. The rapid 
convergence of the variational wave functions with the 
number of Gaussians increasing from N — 1 to 5 is im- 
pressive as can be seen in the insets in Figs. [2] and [3] 
Note that far from the bifurcation the wave function of 



TABLE I. Dependence of the critical scattering length oJc 
of the tangent bifurcation and the mean field energy and the 
chemical potential of the ground state (g) and the excited 
state (e) at scattering length a sc — — 1 on the number N of 
coupled Gaussian functions. 
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-1.178097 
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-1.304592 
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-1.032780 


-0.140151 


-0.584799 


-0.136826 
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-1.025527 


-0.140637 


-0.595154 


-0.138380 


-0.862562 


4 


-1.025167 


-0.140655 


-0.595659 


-0.138449 


-0.860862 
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-1.025149 


-0.140656 


-0.595682 


-0.138452 


-0.860767 


oo 


-1.025147 


-0.140657 


-0.595685 


-0.138453 


-0.860757 
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FIG. 2. Wave function ip{r) of (a) the stable ground state, 
(b) the excited state, at scattering length a sc — — 1 close 
to the exact numerical bifurcation. The calculations with 
coupled Gaussians converge rapidly to the exact numerical 
wave function (triangles). 



the stable state in Fig. [3la) and the wave function of the 
excited state in Fig. [3^b) differ greatly, while for a sc = —1 
both wave functions in Fig. [2|a) and (b) are similar. At 
the tangent bifurcation, both merging states, the stable 
ground state and the excited state are described by the 
same wave function. The reason is that due to the nonlin- 
earity of the GPE, the solution at the tangent bifurcation 
has the properties of an "exceptional point" [5] . 

Both Figs. [5] and [3] apparently show that the form of 
the exact numerical wave function differs from the simple 
N = 1 Gaussian form. The inclusion of five Gaussians, 
however, achieves excellent results. Therefore we can as- 
sume the variational ansatz using five coupled Gaussian 
functions to be converged with sufficient accuracy to cal- 
culate all major properties of the system. 



FIG. 3. Same as Fig. [2]but at scattering length 



-0.6. 




FIG. 4. The constituents of the five-Gaussian ansatz (No. 1, 
. . ., 5) plotted separately. The superposition of the five Gaus- 
sians (coupled), excellently agrees with the exact numerical 
solution (exact) at scattering length a sc — — 1. 



B. Stability of stationary states 



To visualize the ansatz, the five Gaussian functions as 
constituents of the N = 5 solution at scattering length 
a sc = —1 in Fig. [5] are drawn in Fig. |4J The two nearly 
identical curves at the top are the variational solution 
given as the sum of the functions below, and, for com- 
parison, the exact numerical solution (triangles). 



For the spherically symmetric monopolar condensate 
we restrict the discussion of the stability of states to fluc- 
tuations of the wave function in the radial coordinate r. 
In numerically exact calculations, the linearization of the 
GPE with the Frechet derivative [_>] leads to two coupled 
Bogoliubov equations for the real and imaginary parts of 



the wave function ^ Re (r, t) and ip lm (r,t), 
d 
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where ip±(r) is the numerically exact stationary ground 
or excited state, respectively. The method for solving 
those equations with the ansatz for the perturbations 
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where A is one of the exact stability eigenvalues, is elab- 
orated in [6] . The exact stability eigenvalues are used for 
comparisons with the eigenvalues of the Jacobian 
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which is a (4iV x 4.ZV)-dimensional non-symmetric real 
matrix obtained by linearization of the dynamical equa- 
tions of motion Q for the variational parameters pQ. 
The eigenvalues of the Bogoliubov equations ( 12 ) with 
the ansatz ( fl3| and of the Jacobian J in Eq. ( 14 1 always 
occur in pairs ±A with opposite sign. 
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FIG. 5. First pair of eigenvalues as a function of the scat- 
tering length, for an increasing number of coupled Gaussians 
and for the exact numerical calculation, (a) The two low- 
est eigenvalues of the stable ground state are purely imagi- 
nary, (b) The unstable excited state with a pair of purely 
real eigenvalues. Vanishing real or imaginary parts are not 
shown. The insets illustrate the rapid convergence of the vari- 
ationally computed eigenvalues to the exact solutions of the 
Bogoliubov equations. 



1. First pair of stability eigenvalues 

For the ansatz with a single Gaussian there is, after 
exploiting the normalization condition, only one varia- 
tional parameter a — a 1 in Eq. pi for the (complex) 
width of the Gaussian function, and the eigenvalues of 
the Jacobian can be calculated analytically [6l [18] . For 
the ground state the two eigenvalues 



\s _ 



±- 



16i 



1 + 



3tt 



9tt 



1 



8a 3c 
3tt 



1 



(15) 



are purely imaginary for all scattering lengths above the 
critical value of the tangent bifurcation. For the excited 
state there are two purely real eigenvalues 
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For N > 2 coupled Gaussian functions the Jacobian J is 
diagonalized numerically. The stability eigenvalues ob- 
tained with a single Gaussian function (Eqs. (15 1 and 



(16)), the corresponding pair of eigenvalues of the Ja- 
cobian J for TV = 2 to N = 5 coupled Gaussians, and 
the corresponding pair of the exact eigenvalues are pre- 
sented in Fig. [5] Similar to the calculation of the mean 
field energy and chemical potential, the pair of eigenval- 
ues merges and vanishes at a tangent bifurcation. As we 
include more Gaussian functions, the critical scattering 
length shifts to higher values and converges to the bifur- 
cation of the exact numerical solution as was expected 
from the behavior of the energies. Fig. p)la) shows the 
first pair of eigenvalues of the stable ground state. They 
are purely imaginary. Since this is true for all stability 
eigenvalues of the ground state (see below), the branch 
of the ground state is stable. Fig. [5^b) shows the first 
eigenvalue pair of the excited state. These eigenvalues are 
purely real, and thus the excited state is unstable. The 
tangent bifurcation is clearly exhibited for both states, as 
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FIG. 6. (a) Second pair of eigenvalues as a function of 
the scattering length, for varying number of coupled Gaus- 
sians and for the exact numerical calculation for the stable 
ground state and the unstable excited state. Numerically ex- 
act eigenvalues (triangles) are only available for the stable 
ground state. All eigenvalues are purely imaginary, vanishing 
real parts are not shown, (b) Same as (a) but for the third 
pair of eigenvalues. 



each pair of eigenvalues ±A merges at zero and vanishes 
at the critical scattering length. 

It is important to note that similar to the convergence 
properties of the mean field energy and the chemical po- 
tential, very few coupled Gaussians are already sufficient 
to achieve excellent results for the lowest stability eigen- 
values. Results obtained with N > 3 coupled Gaussians 
can not be distinguished in Fig. [5] from the numerically 
exact values. 



2. Additional pairs of stability eigenvalues 

In contrast to the calculation with a single Gaussian, 
which can provide only one pair of eigenvalues, additional 
eigenvalues are accessible when using coupled Gaussian 
wave functions. We compare them with the exact stabil- 
ity eigenvalues in Fig. [6] The additional eigenvalues are 
increasingly difficult to obtain in the numerically exact 
computation. In the following, we refer to the eigen- 
values with the second lowest absolute value simply as 
"second (pair of)" eigenvalues, etc. Numerically exact 



data is available for the lowest three pairs of eigenvalues 
of the stable ground state, and therefore here we com- 
pare the variational solution only with the lowest three 
eigenvalues of the stable solution. 

Figure l(3ta) shows the second pair of eigenvalues for the 
stable ground state and the unstable excited state. The 
second pairs of eigenvalues are all purely imaginary, and 
converge with increasing number of coupled Gaussians to 
the numerically exact eigenvalues. 

The third pair of eigenvalues presented in Fig. pTb) is 
qualitatively similar to the second pair. For any num- 
ber of coupled Gaussians, they are purely imaginary. 
Figs. [6ja) and (b) suggest that, as the number of the 
eigenvalue rises, the convergence progresses more slowly. 
However, using five Gaussians in Fig.^STb) even the third 
pair of the variational eigenvalues shows no apparent de- 
viation from the numerically exact result. 



3. Variations of the ground state wave function 

We investigated the stability of the stationary states by 
linearizing the dynamical equations in the vicinity of the 
fixed points. For the stable ground state there are only 
purely imaginary eigenvalues, for the unstable excited 
state we found one pair of real eigenvalues. Additionally 
to the analysis of the eigenvalues \ of the Jacobian of 
the linearized dynamical equations, we can evaluate the 
respective eigenvectors, which provide the form of the 
wave function's fluctuations. 

We focus on variations of the Gaussian parameters cor- 
responding to the eigenvector % and calculate the first or- 
der power series of the wave function ip(r, t) at the fixed 
point (FP) p, 
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dy)i(r,t) = y j \\r oa^ — r oa i +107,^ — 07* 
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x/rMe 
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For a scattering length of a sc = —0.8, Fig. [7] shows the 
variation of the ground state wave function for the eigen- 
vectors corresponding to the first (a), second (b) and 
third (c) pair of eigenvalues ±A for the variational calcu- 
lation, as well as for the numerically exact calculation. 

In Fig. [7[a) the variation of the wave function Sip con- 
verges rapidly with increasing number of Gaussian func- 
tions to the numerical solution. For the second and third 
pair of eigenvalues in Fig.^b) and (c), the variation us- 
ing five Gaussians is almost identical to the numerical 
variation. For the higher pairs of eigenvalues, the solu- 
tion obviously converges slower. However, we note that it 
also becomes more and more difficult to obtain the exact 
solutions. The spatial extension of the fluctuations Sip 
exceeds the elongation of the wave function ip (cf. Figs. [2] 
and p5| by far and it becomes hard to achieve converged 
fluctuations. The differences between the numerical and 
the N = 5 solution in Fig. Wic) can already be explained 
by the quality of the numerical approach. 



(a) 4.: 



N = 2 
N = 3 
N=4 
N = 5 

numerical 




2 4 6 



10 12 14 



FIG. 7. Deviations of the ground state wave function ac- 
cording to the eigenvectors of the (a) first, (b) second, and 
(c) third pair of eigenvalues ±A of the stability analysis at 
scattering length a BC — —0.8. 



III. DIPOLAR CONDENSATES 

The stationary GPE for a dipolar BEC with a har- 
monic trap, the short-range s-wave scattering term and 
the long-range dipolar interaction reads 



- A + -fix 2 + rfy 2 + j 2 z 2 + 8ttA--- \ip(r) 
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In Eq. (TT8 ) the "natural units" for this system intro- 
duced in |14| have been used, which are h for action, 



rrid = 2m for mass, ad = mdfJ.ol^ 2 /(^fi 2 ) for length, 
Ed = h 2 /(mda d ) for energy, and uod = Ed/fr for fre- 
quency. The angle between the external magnetic field 
in z-direction and the vector r — r' is denoted 9, and N 
is the particle number. As in 14 we scale r = Nf, 
if) = N~ 3 > 2 ip, insert the newly defined quantities in 
Eq. (jT8l), redefine % = N 2 ji, p, = N 2 fi, and after- 



wards omit the tilde once again. With the replacement 
fi — > i(d/di) we finally obtain the time-dependent GPE 
for a dipolar BEC in particle number scaled dimension- 
less units, 



A + -fix 2 + 7,V + 7 fz 2 + 87ra sc |^(r)| z 
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with the trap frequencies J x ,y,z = N 2 uj xyz /(2u)d) and 
the s-wave scattering length a sc — aj ad- For dipolar con- 
densates it is possible to solve the GPE fully numerically 
on a two- or three-dimensional lattice [TSJ HI] • Ronen et 
al. [TT] and Dutta et al. [E5] have shown that in certain re- 
gions of the parameter space dipolar condensates assume 
a non-Gaussian, biconcave, "blood-cell-like" shape. In 
this paper we want to apply the variational method with 
coupled Gaussian functions introduced in the preceding 
paper [T]. We will show that the variational technique 
is a full-fledged alternative to the numerical simulations 
on grids, and additionally uncovers unstable stationary 
solutions not accessible in previous full-numerical evalu- 
ations. 

The frequency and symmetry of the magnetic trap 
strongly influences the physical behavior of dipolar Bosc- 
Einstcin condensates. In the following we will analyze 
one distinct trap symmetry in detail, where the conden- 
sate has a blood-cell-shaped form. The ansatz with cou- 
pled Gaussians has proved to be capable to modify the 
simple Gaussian form of the wave function for monopolar 



condensates (see Sec. II A). The biconcave shape, how- 



ever, where the maximum density is no longer located at 
the origin is even a stronger challenge for the variational 
approach. We investigate the dipolar condensate for an 
axially symmetric trap with trap frequencies 

lx = ly = le = 3600, lz = 25200, 

which is equivalent to the frequency ratio and 
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and corresponds to a value of D = y/j^/2 = 30 in [llj . 
The trapping frequency in the z-direction parallel to the 
orientation of the dipoles is seven times larger than in 
the plane perpendicular to that direction, and for some 
parameters a sc the ground state of the condensate has 
a biconcave, blood-cell-shaped form [TT] . In contrast to 
monopolar condensates where the inclusion of additional 
Gaussian functions provides an improved numerical ac- 
curacy of the results, dipolar condensates offer a wealth 



of new phenomena with increasing number of coupled 
Gaussian functions as will be shown below. 



A. Variational calculations with one and two 
Gaussian functions 

As ansatz for the variational calculations we use the 
wave function 



N N 

v,(r-,i) = E ei(a " x2+a ^ 2+a ^ 2+7fc) ^E-' 



(20) 



fc=i 



fc=i 



where N is the number of coupled Gaussians. For an 
axisymmetric trap the stationary solutions are also sym- 
metric, i.e., a^ = aJl = a' k . Nevertheless, all stabil- 
ity properties have been computed with the fully three- 
dimensional ansatz. The case of a single Gaussian func- 
tion (N = 1) has been discussed in [14] . In this section 
we demonstrate that results especially for the mean field 
energy and chemical potential are already substantially 
improved with the use of only N = 2 coupled Gaussians. 
Results of variational calculations with a single Gaussian 
function and numerical simulations on grids are employed 
for comparison and discussion. 
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1. Stationary solutions 



Using the variational ansatz (20 1 stationary solutions 



of the Gross-Pitaevskii equation for dipolar BEC are ob- 
tained by a numerical root search for the fixed points of 
the dynamical equations for the variational parameters 
as described in detail in the preceding paper pQ. 

Figure p^ presents in (a) the mean field energy and in 
(b) the chemical potential for the two stationary solutions 
obtained with a single Gaussian function and with two 
coupled Gaussians. For the ground state results of a nu- 
merical lattice calculation are also marked in Fig. [8] The 
numerical simulation was performed on a lattice with a 
grid size of 128 x 512 points using fast-Fourier techniques 
and imaginary time evolution of an initial wave function. 

The variational calculations with one Gaussian (N = 
1) show the following behavior. For scattering lengths 
below a^' var = -0.0378917 there is no stable conden- 
sate. Similar as in monopolar condensates two solutions 
are born at the critical scattering length in a tangent 
bifurcation, the stable ground state (vl) and an unsta- 
ble excited state (v2). For a detailed stability analysis 
see Sec. MI A 21 below. The unstable branch vanishes at 
scattering length a sc = 1/6. 

The variational ansatz with N = 1 is limited to the 
Gaussian shape of the wave function with two width pa- 
rameters a e and a z , and thus the values obtained for the 
mean field energy and chemical potential are not very ac- 
curate. However, the results are substantially improved 
when using a variational ansatz with two coupled Gaus- 
sians. This can be seen in Fig.|8]especially for the ground 



FIG. 8. (a) Mean field energy and (b) chemical potential as a 
function of the scattering length, variational solution (v) using 
two Gaussians compared with an exact full-numerical lattice 
calculation (nl). For two Gaussians (TV = 2), the ground 
state (vl) and the excited state (v2) emerge in a tangent 
bifurcation at a sc — —0.034202. The unstable branch of the 
chemical potential in (b) has a maximum at /i max ~ 170000 
(not shown) and then vanishes as a nearly vertical line at 
a BC — 1/6. 



state when comparing the N — 2 variational computa- 
tion with the lattice computation (nl) marked by the 
triangles in Fig. [8] 

In the full-numerical grid calculations only the ground 
state can be obtained. Starting with positive scatter- 
ing lengths and decreasing a sc , the numerical grid calcu- 
lations provide a ground state down to a critical point 
a cr,num = _q.008. Note that the ground state of the 
solution using two coupled Gaussian wave functions is 
nearly indistinguishable from the numerical lattice calcu- 
lation in Figs. (8la) and (b). An important advantage of 
the variational method is that it can provide both stable 
and unstable states. As will be shown below, the stabil- 
ity properties of the variational solutions can clarify the 
mechanism of how the condensate turns unstable. 

Regarding the wave functions, one single Gaussian can 
evidently not adequately represent a blood-cell-shaped 
condensate. Two Gaussians not only significantly in- 
crease the accuracy of the mean field energy, but also 
greatly improve the form of the wave function. Even the 
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FIG. 9. Eigenvalues of the Jacobian for the variational solu- 
tion with one Gaussian as a function of the scattering length 
a sc . (a) The ground state is stable, all eigenvalues are purely 
imaginary, (b) the excited state is unstable since there are 
real parts of eigenvalues, emerging in a tangent bifurcation 
at aJJ = —0.0378917. Eigenvalues which do not reach A = 
at Usl match with the corresponding eigenvalues of the stable 
and unstable state, respectively. 



FIG. 10. Same as Fig. [9] but for the variational calculation 
with two coupled Gaussians, (a) the stable ground state and 
(b) the unstable excited state. Additional eigenvalues are re- 
vealed. The tangent bifurcation is shifted to aj£ = —0.034202. 
Eigenvalues which do not reach A = at a^l match with the 
corresponding eigenvalues of the stable and unstable state, 
respectively. 



biconcave shape of the dipolar condensate is qualitatively 
visible in the variational solution with two-Gaussians, 
however, the result is not fully converged. Exact wave 
functions will be compared in Sec. |IIIB] with variational 
results obtained with more than two coupled Gaussians. 



2. Stability analysis 

To perform a stability analysis with numerical lattice 
calculations the Bogoliubov-de Gennes equations have to 
be solved [12 [2U] . Here we restrict our discussion to the 
stability analysis of the variational solutions, which is in- 
structive considering nonlinear dynamics and bifurcation 
theory. 

We follow the procedure outlined in pQ and start with 
the stationary solutions of the GPE calculated using one 
and two Gaussians. These solutions are fixed points of 
the dynamical equations for the Gaussian parameters. 
We then linearize these dynamical equations in the vicin- 
ity of the fixed point and calculate the eigenvalues of the 
Jacobian (see pQ). 

The eigenvalues of the ground state obtained with one 



Gaussian wave function in Fig. K3[a) are purely imagi- 
nary. Therefore this state is stable. If we perturb the 
variational parameters of the fixed point solution, the 
quasi-periodic motion is confined to the vicinity of the 
fixed point. Figure|9{b) shows the characteristic eigenval- 
ues of the unstable excited state. Contrary to the stable 
state, there are eigenvalues with non-vanishing real parts 
Re A. Perturbations in the direction of the corresponding 
eigenvector lead to an exponential growth of the pertur- 
bation. Therefore this state is unstable. Both branches 
exist for scattering lengths down to af c — —0.0378917, 
where they merge in a tangent bifurcation (see Fig. M) . 
The tangent bifurcation is apparent in the eigenvalues in 
Figs. [9|a) and (b). 

Some eigenvalues of the Jacobian using a wave func- 
tion of two coupled Gaussians in Fig. [10] qualitatively 
agree with the one-Gaussian calculation. There is one 
stable ground state with purely imaginary eigenvalues in 
Fig. [l0|a), and one unstable excited state in Fig. |l0|(b). 
However, the two-Gaussian calculation yields additional 
eigenvalues which are not available within the limited pa- 
rameter space of the simple one-Gaussian calculation. 

The unstable branches of both calculations in 



10 



Figs. [9Vb) and[l0|b), exhibit that for the given parame- 
ters of the trap in dipolar condensates, the stability sce- 
nario is quite more complex than for monopolar conden- 
sates. As we can see, there is not only one single pair of 
imaginary eigenvalues of the stable solution approaching 
zero and merging with one pair of real unstable eigen- 
values (denoted 1 in Fig. [9lb)) from the excited state. 
For monopolar condensates this real pair of eigenvalues 
of the unstable solution remains the only one for increas- 



ing scattering lengths (see Sec. II B). By contrast, here, 
a second pair of eigenvalues (denoted 2 in Fig. |9jb)) in- 
dicating instability additionally forms as we follow the 
excited state from the bifurcation point to positive scat- 
tering lengths. The eigenvectors that correspond to this 
pair of real eigenvalues also show an interesting behav- 
ior: The stability analysis is performed in three dimen- 
sions although the solutions of the GPE are axisymmet- 
ric. Therefore in the linearization around the fixed point, 
perturbations in x and y direction are calculated indepen- 
dently. The eigenvectors corresponding to the additional 
unstable eigenvalues are not symmetrical in x and y, and 
thus break the axial symmetry of the fixed point solution. 
The unstable branch of the calculation with two Gaus- 
sians already shows multiple pairs of unstable real eigen- 
values. This suggests that the dynamics of dipolar con- 
densates in this parameter region is very complex and 
can be described better by including even more coupled 
Gaussians in the variational approach. Indeed, the fur- 
ther increase of the number of variational parameters re- 
veals new physical properties and phenomena also of the 
ground state of the biconcave dipolar condensate. 



B. Converged variational calculations with up to 
six coupled Gaussian functions 
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FIG. 11. (a) Convergence of the mean field energy (to the 
solid line) with increasing number of coupled Gaussian wave 
functions for a sc — 0. The mean field energy for four cou- 
pled functions lies already energetically lower than the nu- 
merical value of the lattice computation (dashed line), (b) 
The converged wave function ip as a function of the trans- 
verse coordinate g at different z coordinates. The variational 
solution and the numerical lattice solution are denoted v and 
n, respectively. 



Since the inclusion of two Gaussians already substan- 
tially improves the mean field energy, the coupling of 
more functions only results in a minor correction in the 
value of the mean field energy, which would not be ap- 
parent in, e.g., Fig. [8J We therefore present in Fig. 11 'a) 
the convergence of the mean field energy at one selected 
scattering length, viz. a sc = 0, for N = 2 to N = 6 
coupled functions. For other scattering lengths the con- 
vergence behavior is similar. This example shows that 
as few as four coupled Gaussian functions result in a 
mean field energy which lies below the numerical solu- 
tion of the lattice calculation (dashed line) . For five and 
six Gaussians, the variational solution converges to a dis- 
tinct value (solid line). Note that the simplest variational 
solution with one Gaussian function is not included in the 
figure because the mean field energy of E m { — 60361 lies 
far outside the vertical energy scale. 

In Fig. fTTVb) for the same scattering length (a sc = 0) 
the converged wave function is shown at different z coor- 
dinates. The wave function of the variational calculation 
is practically identical to the wave function of the nu- 
merical lattice calculation. Both wave functions show 



the characteristic biconcave shape of the condensate. 

In this section, we discuss properties of the solutions 
obtained with five and six coupled Gaussians, which are 
qualitatively identical and quantitatively indistinguish- 
able in the figures presented. Therefore we will omit the 
detailed label and refer to the converged variational solu- 
tion simply as "variational solution" (v). While the use 



of one and two coupled Gaussians in Sec. Ill A results 
in two branches, one stable, one unstable, emerging in a 
tangent bifurcation, this bifurcation scenario has to be 
revised with the converged variational solution. 

Figure |12[ a) shows an overview of the mean field en- 
ergy for the variational solution. There are two impor- 
tant intervals of the scattering length a sc , showing differ- 
ent characteristics of the variational solution with cou- 
pled Gaussian functions. These intervals of the scattering 
length are marked in Fig. 12 ^a). The different line styles 
in Fig. [12^ a) indicate the stability of the solutions antici- 
pating the results of the stability analysis. The numerical 
solution via lattice calculation and imaginary time evo- 
lution obtains only the ground state. Note that the nu- 
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FIG. 12. (a) Overview of the mean field energy for the bi- 
furcation scenario revealed in the variational calculation. For 
comparison some values of the numerical lattice solution for 
the ground state are presented (see text for more details). The 
stability eigenvalues in the regions around the critical scatter- 
ing lengths a£'* = -0.01224 and a c / c ,p = -0.00359 marked by 
the black frames in (a) are shown in (b) and (c), respectively. 
In (b) I and II denote the two merging branches. Only the 
eigenvalues involved in the bifurcation are shown (for analysis 
and interpretation see text). Eigenvalues of the Jacobian as a 
function of the scattering length in (c) indicate the stability 
change of the ground state in a pitchfork bifurcation. The un- 
stable state for a sc < asc' p — —0.00359 turns into the stable 
ground state for a sc > a^' p in a pitchfork bifurcation. Sub- 
figure (c) only shows the lowest eigenvalues, those involved in 
the stability change. 



merical simulation was carried out on a two-dimensional 
axisymmetric grid, and thus the imaginary time evolu- 
tion can provide unstable ground states if the instability 
is rotational, i.e., resulting in an angular collapse of the 
condensate [T5] . 

The variational solution is able to obtain both, the sta- 
ble ground state and stationary excited states. There 
are variational results down to a critical point a* r = 
—0.01224. To analyze the stability of the solution, we 
present in Fig. [l2J[b) a stability analysis of the linearized 
dynamical equations in the interval —0.0123 < a sc < 
—0.0088 of the scattering length [frame marked (b) in 
Fig. Ilia)]. 



Figure fT2jb) shows at the center the typical scenario 
of eigenvalues of two branches merging in a tangent bi- 
furcation at a* r = —0.01224. A pair of purely imaginary 
eigenvalues of the branch denoted I merges at a critical 
point with a pair of real eigenvalues of branch II. Respec- 
tive vanishing real or imaginary parts are not shown. 

In addition to this tangent bifurcation there is a direc- 
tion in Gaussian parameter space, in which both branches 
show unstable, purely real eigenvalues (see Fig.[l2^b) top 
and bottom). Therefore both branches involved in the 
tangent bifurcation are born unstable. 

The previous scenario which resulted from the calcula- 
tion with one or two coupled Gaussians in Sec. |IIIA| must 
now be revised. In the converged variational ansatz, the 
tangent bifurcation is on top of an unstable direction. 
The important conclusion is that there is no stable con- 
densate in this region of the scattering length. 

Where does the variational condensate turn stable? To 
pursue this question we increase the scattering length to 
a sc s=s -0.00359 [frame marked (c) in Fig. [l2^a)]. The 
corresponding stability analysis shows a stability change 
for the lowest eigenvalues of the ground state, which are 



plotted in Fig. 12 c) 



For scattering lengths a sc < af c ' p = —0.00359 the 
branch is unstable, indicated by the pair of real eigen- 
values in Fig. fl2tc). At this bifurcation point the real 
eigenvalues vanish and for a sc > a^ ,p = —0.00359 a sta- 
ble ground state forms, indicated by a pair of imaginary 
eigenvalues. Figure 12 [c) shows only the respective low- 



est pair of eigenvalues which is involved in the stability 
change, all other eigenvalues are purely imaginary, and 
arc omitted for the sake of clarity of the figure. 

The stability change of the ground state in Fig. [l2^a) 
takes place in a pitchfork bifurcation. From left to right, 
one unstable branch turns stable in the bifurcation. 

In general, three branches are involved in a pitchfork 
bifurcation [21] . If the ground state of the dipolar con- 
densate changes stability in a pitchfork bifurcation, two 
stable states should appear as stationary states in the 
variational calculation for the dipolar BEC as well. How- 
ever, for the following reasons we are probably unable to 
observe these additional stable states directly. There are 
AN complex Gaussian variational parameters, i.e. 48 real 
parameters for six coupled functions. The pitchfork bi- 
furcation and the stability change take place in one direc- 
tion characterized by the eigenvectors of the eigenvalues 
shown in Fig. |12|^c). Since an increasing number of vari- 
ational parameters leads to a more and more complex 
parameter space with increasingly complex interactions 
between the degrees of freedom, it is well possible that 
the two stable states are limited to an extremely tiny 
vicinity of the bifurcation point. 

For scattering lengths greater than the bifurcation, the 
ground state is stable. At the bifurcation the system 
also changes from regular dynamics to chaotic dynam- 
ics, where for scattering lengths below the bifurcation 
point both additional stable branches may undergo sev- 
eral bifurcations themselves, immediately turning them 
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FIG. 13. Contour plot of the mean field energy (a) for 
«sc = —0.0036 closely below, and (b) for a sc = —0.00358 
closely above the pitchfork bifurcati on. The eigenvectors cor- 
responding to the eigenvalues in Fig. 12 c) linearize the vicin- 
ity of the fixed point (Si and $2, arbitrary units). 



unstable. Therefore it is not possible to obtain the sta- 
ble branches directly via search for fixed points of the 
dynamical equations. However, it is possible to catch a 
glimpse of those branches, if we consider the linearized 
surroundings of the stability changing state very close to 
the bifurcation. 

The eigenvectors corresponding to the lowest eigenval- 
ues (that show the stability change) linearize the vicinity 
of the stationary state (61,62)- Figure 13 shows the con- 



tour plot of the mean field energy of this linearization in 
arbitrary units for two scattering lengths (a) very close 
below and (b) very close above the bifurcation. Above 
the bifurcation Fig. 13 ]b) shows one elliptic stable state, 
the stable ground state. Below the bifurcation, Fig. 13 'a) 
shows the unstable fixed point at the center. Besides the 
unstable hyperbolic fixed point, there are two stable el- 
liptic points at (±0.5,^0.5) in this vicinity linearized by 
the eigenvectors. Nevertheless, Fig. [13] is limited to the 
two-dimensional plane spanned by two eigenvectors. If 
all directions of the eigenvectors of all eigenvalues are 
considered, those two states are only stable in a very 
small interval a" c ,p — e < a sc < af c - p below the bifur- 
cation point, which makes them numerically impossible 
to find. Due to the numerically small value of e, how- 
ever, the classification of the condensate as unstable for 
scattering lengths a sc < a p r — —0.0359 remains true in 



physical terms. 

If we further investigate the two eigenvectors that cor- 
respond to the pair of eigenvalues in which the stability 
change occurs, we see, that the axial symmetry is no 
longer present. For the present trap symmetry and fre- 
quencies (j x — j y = 3600 and j z — 25200) and the 
ansatz 



A? 



V^,z)=£ ei(a>2+ ^ 2+7fc) ' 



(21) 



fc=i 



all fixed points had been axisymmetric. The stability 
analysis, however, is done without any assumptions con- 
sidering symmetry, allowing variations in both, 6a x and 
6a 1 ! separately. Therefore, oscillations in directions which 
break the axial symmetry are allowed. 

The characteristic eigenvectors can be considered as 
deviations of the wave function Sip (see [1]). If the eigen- 
vector is no longer axially symmetric, i.e., 5a^ — — Sa^ 
for all k, the perturbation leads to an asymmetric oscil- 
lation or collapse of the condensate. Indeed, we find this 
kind of eigenvectors for the lowest eigenvalues of the Ja- 
cobian for the variational solution of the ground state in 
Fig. [l2J[c) . This behavior of the eigenvectors of the Ja- 
cobian is an indication of the so called "angular collapse 
of dipolar BEC" associated with the biconcave shape of 
the condensate 13 . 

This angular collapse can be observed in a time evolu- 
tion of the condensate. We prepare the stationary wave 
function for a scattering length closely below the bifur- 
cation, and add a deviation 



ip(z) = ip(z 



FP\ 



■ 5ip(6zi) 



(22) 



in the direction of the eigenvector i whose corresponding 
eigenvalue is involved in the stability change. Figure 14 
shows a time evolution of the particle density |^| 2 ob- 
tained by numerical integration of the dynamical equa- 
tions. The time evolution of the unstable excited state 
clearly reveals an angular collapse of the condensate, the 
particle density concentrates on two non-axially symmet- 



ric regions as shown in Fig. 14 The Gaussian ansatz ( 20 ) 



implies a collapse of the condensate with parity with re- 
spect to the x and y axis, but it may be possible in fu- 
ture work to modify the Gaussian functions to include 
any form of collapse without symmetry restrictions. A 
modified ansatz may even allow for an angular collapse 
of the condensate with three density peaks as reported 
by Wilson et al. [13]. 



IV. CONCLUSION 

We have applied the method of coupled Gaussian wave 
packets to Bose-Einstein condensates with two different 
types of long-range interaction, viz. an attractive gravity- 
like 1/r interaction and a dipole-dipole interaction. The 
mean field energy and chemical potential have been ob- 
tained as fixed points of dynamical equations for the set 
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FIG. 14. Time evolution of the perturbed particle density 
\ip\ 2 of the unstable stationary solution as a function of x and 
y for z = for (a) t = 0.0001, (b) t = 0.005, and (c) t = 0.006. 



of variational parameters. As an alternative to solving 
the Bogoliubov equations the stability properties of the 
condensates have been determined by applying methods 
of nonlinear dynamics to the linearized equations of mo- 
tion. 

For monopolar condensates we have shown that the ad- 
ditional variational parameters of the coupled Gaussian 
ansatz greatly improve the accuracy of the variational so- 



lution in comparison to the established single Gaussian 
ansatz. With three coupled Gaussian functions in the 
trial wave function, the numerical mean field energy is 
already reproduced with an accuracy of more than four 
digits. The solution with five Gaussians proves to be fully 
converged to the solution of the direct numerical integra- 
tion of the GPE. Furthermore, the stability properties 
and the bifurcation of the numerical solution are excel- 
lently reproduced by the coupled Gaussian ansatz. The 
variational method also provides easy access to higher 
stability eigenvalues, which numerically are hard to ob- 
tain. For monopolar condensates, the method of coupled 
Gaussian functions is an excellent and fully valid alter- 
native to the direct numerical integration of the GPE. 

For dipolar condensates we have described the new phe- 
nomena revealed by variational solutions with an increas- 
ing number of coupled Gaussians. The variational ansatz 
with multiple coupled Gaussian functions turns out to be 
a full-fledged alternative to numerical lattice calculations 
for condensates with dipolar interaction. With the use of 
as little as five to six Gaussian functions, the variational 
solution can be considered to be fully converged. In con- 
trast to lattice calculations via imaginary time evolution, 
the variational ansatz also obtains excited states. Thus 
the method of coupled Gaussian functions gives access 
to yet unexplored regions of the space of solutions of the 
GPE, and we have been able to clarify the theoretical 
nature of the collapse mechanism: The ground state of 
the condensate turns unstable in a pitchfork bifurcation 
before it finally vanishes in a tangent bifurcation. The 
stability analysis indicates a further feature of the col- 
lapse mechanism: The condensate breaks the cylindrical 
symmetry on the verge of collapse, indicating an angular 
decay of the condensate. 

The convergence of the variational method with Gaus- 
sians has proved to be very fast even close to the crit- 
ical scattering length, at which the collapse of the con- 
densate sets in. In future work it will be interesting to 
monitor the convergence of the ansatz in the Thomas- 
Fermi regime, where exact polynomial solutions of the 
wave functions are found -2, 22, .23 J. It will also be desir- 
able to investigate in more detail how the numerical sta- 
bility analysis via the Bogoliubov-de Gennes equations 
is related to the variational stability analysis. Further- 
more it will be possible to investigate real time dynamics 
of the decay of dipolar BEC and angular rotons with a 
modified coupled Gaussian ansatz. The coupled Gaus- 
sian ansatz proved to be a fast and accurate alternative 
to full-numerical calculations. There are several interest- 
ing systems where this method can be applied in future 
work, e.g., monopolar BEC with vortices or stacks of 
dipolar condensates. 
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